with(combinat):
with(ListTools):
with(Iterator):

print(`To see the procedures in this package, type "Help()"`):
print(`For explanation of each procedure, type "WhatIs(procedure_name)"`):

Help:=proc(): print(`RanDeck(n,k,r)`): 
print(`WARm(m,Q,r)`, `WARmax(m,Q,r)`, `WARmin(m,Q,r)`):
print(`RepPlayed(Q,r,b,N)`, `RepAllAsc(Q,r,b,N)`,`RepAllDesc(Q,r,b,N)`, `RepAsc(Q,r,b,N)`, `RepDesc(Q,r,b,N)`):
print(`GenOneMove(P,c,m,WAR,rep)`, `GenOneGameK(n,k,r,m, WAR,rep,K)`, `GenManyGamesK(n,k,r,m,WAR,rep,K,w,t,M)`, `GenManyGamesKEDIT(n,k,r,m,WAR,rep,K,w,t,M)`, `
FinishGame(H,r,m, WAR,rep)`):
print(`SimpPerDecks(n,rep)`, `SimpGameLengths(n,rep)`, `SimpGameInfo(n,rep)`):
end:

WhatIs:=proc(procedure) :
  if procedure=RanDeck then
    print(`inputs positive integers n,k,r`):
    print(`Simulates dealing out as many cards from a shuffled deck containing k copies of the cards {1,...,n} to r players such that every player has the same number of cards.`):
    print(`Outputs a list of length r, where the i-th entry is a list of length floor(n*k/r) that represents the i-th players deck.`):
 
   elif procedure=WARm then
     print(`inputs a non-negative integer m, a list Q, and a positive integer r`):
     print(`m represents the number of cards placed face down in a WAR`, 
		`Q is each players deck at the start of the WAR (not including the cards that triggered the war)`, 
		`r is the number of players at the start of the game (should have nops(Q)=r, where players who have lost have the deck [])`):
     print(`Simulates a WAR where m cards must be placed face down NO MATTER WHAT. Any player that does not have enough cards automatically loses.`):
     print(`If no one can complete the WAR it outputs {0}, indicating that the game will have no winners.`):
     print(`Otherwise, it outputs a list [tempQ,b], where tempQ is a list of each players deck AFTER the war and b is a list of the cards that were either placed face down or played during the WAR`):

   elif procedure=WARmin then
     print(`inputs a non-negative integer m, a list Q, and a positive integer r`):
     print(`m represents the number of cards placed face down in a WAR`, 
		`Q is each players deck at the start of the WAR (not including the cards that triggered the war)`, 
		`r is the number of players at the start of the game (should have nops(Q)=r, where players who have lost have the deck [])`):
     print(`Simulates a WAR where m cards must be placed face down UNLESS any player still in the game does not have enough cards to do so.`, 
		`If the smallest hand at the START of the WAR is k<m+1, then every player places k-1 cards face down and plays their next card.`):
     print(`If NO player has any remaining cards it outputs 0, indicating a tie in the war`):
     print(`Otherwise, it outputs a list [tempQ,b], where tempQ is a list of each players deck AFTER the war and b is a list of the cards that were either placed face down or played during the WAR`):

   elif procedure=WARmax then
     print(`inputs a non-negative integer m, a list Q, and a positive integer r`):
     print(`m represents the number of cards placed face down in a WAR`, 
		`Q is each players deck at the start of the WAR (not including the cards that triggered the war)`, 
		`r is the number of players at the start of the game (should have nops(Q)=r, where players who have lost have the deck [])`):
     print(`Simulates a WAR where m cards must be placed face down UNLESS no player has enough cards to do so.`, 
		`If the largest hand at the START of the WAR is k<m+1, then every player places k-1 cards face down and plays their remaining card.`,
		`Any player that does not have enough cards to do so automatically loses.`):
     print(`If NO player has any remaining cards it outputs 0, indicating a tie in the war`):
     print(`Otherwise, it outputs a list [tempQ,b], where tempQ is a list of each players deck AFTER the war and b is a list of the cards that were either placed face down or played during the WAR`):

   elif procedure=RepPlayed then
     print(`inputs a list Q, a positive integer r, a list b, and a positive integer N`):
     print(`Q denotes each players hand after playing all cards in a round`, 
		`r denotes the number of players at the start of the game`, 
		`b is a list of the cards on the table in the order they were played (i.e. the cards that the winner of the round will take)`,
		`Player N is the winner of the round`):
    print(`Outputs a list of each players hand after adding the won cards to the bottom of player N's deck in the order that they were played.`):
    print(`NOTE: order played- player 1 plays a card, then player 2, and so on...`,
		`In a WAR, player 1 places 1 card face down, then player 2,...,then player r, then player 1 places their 2nd card face down,....`):

  elif procedure=RepAllAsc then
     print(`inputs a list Q, a positive integer r, a list b, and a positive integer N`):
     print(`Q denotes each players hand after playing all cards in a round`, 
		`r denotes the number of players at the start of the game`, 
		`b is a list of the cards on the table in the order they were played (i.e. the cards that the winner of the round will take)`,
		`Player N is the winner of the round`):
    print(`Outputs a list of each players hand after sorting ALL the won cards in ascending order and then placing them at the bottom of Player N's deck.`):

  elif procedure=RepAllDesc then
     print(`inputs a list Q, a positive integer r, a list b, and a positive integer N`):
     print(`Q denotes each players hand after playing all cards in a round`, 
		`r denotes the number of players at the start of the game`, 
		`b is a list of the cards on the table in the order they were played (i.e. the cards that the winner of the round will take)`,
		`Player N is the winner of the round`):
    print(`Outputs a list of each players hand after sorting ALL the won cards in descending order and then placing them at the bottom of Player N's deck.`):

  elif procedure=RepAsc then
     print(`inputs a list Q, a positive integer r, a list b, and a positive integer N`):
     print(`Q denotes each players hand after playing all cards in a round`, 
		`r denotes the number of players at the start of the game`, 
		`b is a list of the cards on the table in the order they were played (i.e. the cards that the winner of the round will take)`,
		`Player N is the winner of the round`):
    print(`Outputs a list of each players hand after sorting each subround (i.e. each time all the players placed a card on the table) in ascending order and then placing them at the bottom of Player N's deck.`):
    print(`So, in a WAR, we take the cards each player placed face down first and sort those in ascending order, then take the cards each player played next and sort them in ascending order, and so on.`):

  elif procedure=RepDesc then
     print(`inputs a list Q, a positive integer r, a list b, and a positive integer N`):
     print(`Q denotes each players hand after playing all cards in a round`, 
		`r denotes the number of players at the start of the game`, 
		`b is a list of the cards on the table in the order they were played (i.e. the cards that the winner of the round will take)`,
		`Player N is the winner of the round`):
    print(`Outputs a list of each players hand after sorting each subround (i.e. each time all the players placed a card on the table) in descending order and then placing them at the bottom of Player N's deck.`):
    print(`So, in a WAR, we take the cards each player placed face down first and sort those in descending order, then take the cards each player played next and sort them in descending order, and so on.`):

  elif procedure=GenOneMove then
     print(`inputs a list P, a list c, a non-negative integer m, WAR, and rep`):
     print(`P is each players deck at the start of this turn`,
		`c is a list of the cards that have already been placed on the table during this round (i.e. cards accumulated during unresolved WARs)`,
		`m is the default number of cards to be placed face down during a WAR`,
		`WAR indicates which WAR procedure being used - either WARm, WARmax, or WARmin`,
		`rep is which replacement procedure which replacement procedure is being used - either RepPlayed, RepAllAsc, RepAllDesc, RepAsc, or RepDesc`):
     print(`If this round ends the game with NO winners it outputs {0}`,
		`If this round ends the game with a TIE, then it outputs a set of the players that tied`,
		`Otherwise, it outputs a list Q of each players hand at the end of the round after adding the won cards to the bottom fo the winners deck.`):
		

   elif procedure=GenOneGameK then
     print(`inputs positive integers n,k,r, non-negative integer m, WAR, rep, and a positive integer K`):
     print(`Simulates a game of war that is at most K rounds and has n different cards, k copies of each card, and r players.`,  
		`m is the default number of cards to be placed face down during a WAR`,
		`WAR indicates which WAR procedure being used - either WARm, WARmax, or WARmin`,
		`rep is which replacement procedure which replacement procedure is being used - either RepPlayed, RepAllAsc, RepAllDesc, RepAsc, or RepDesc`,
		`K is the number of rounds we allow before just terminating the game`):
     print(`If the game finishes within K rounds, outputs a list [w,i] where w is the winner (or the set of winners if there is a tie), and i is the number of rounds the game took`):
     print(`If there are no winners, w={0}`):

   elif procedure=GenManyGamesK then
     print(`inputs positive integers n,k,r, non-negative integer m, WAR, rep, a positive integer K, variables w and t, and a positive integer M`):
     print(`Simulates M games of war where each game is at most K rounds and has n different cards, k copies of each card, and r players.`,  
		`m is the default number of cards to be placed face down during a WAR`,
		`WAR indicates which WAR procedure being used - either WARm, WARmax, or WARmin`,
		`rep is which replacement procedure which replacement procedure is being used - either RepPlayed, RepAllAsc, RepAllDesc, RepAsc, or RepDesc`,
		`K is the number of rounds we allow before just terminating the game`):
     print(`outputs a list [W,T,nw,I] where W[i] is the G.f for duration if player i won, T[i] is the G.f for the duration if player i tied with someone else, nw is the G.f for duration if no winners, and I is the number of games that did not end in K moves`):

   elif procedure=GenManyGamesKEDIT then
     print(`inputs positive integers n,k,r, non-negative integer m, WAR, rep, a positive integer K, variables w and t, and a positive integer M`):
     print(`Simulates M games of war where each game is at most K rounds and has n different cards, k copies of each card, and r players.`,  
		`m is the default number of cards to be placed face down during a WAR`,
		`WAR indicates which WAR procedure being used - either WARm, WARmax, or WARmin`,
		`rep is which replacement procedure which replacement procedure is being used - either RepPlayed, RepAllAsc, RepAllDesc, RepAsc, or RepDesc`,
		`K is the number of rounds we allow before just terminating the game`):
     print(`outputs a list [all, win,nw,I] where all is the G.f for duration of all the games that ended, win is the G.f for the duration if 1 winner,  nw is the G.f for duration is no unique winner, and I is the number of games that did not end in K moves`):

elif procedure=FinishGame then
     print(`inputs list H, positive integers r where H has lenght r, non-negative integer m, war procedure WAR, replacement procedure rep`):
     print(`Simulates the rest of the game of war where the i-th entry of H represents the i-th players deck.`,
		`m is the default number of cards to be placed face down during a WAR`,
		`WAR indicates which WAR procedure being used - either WARm, WARmax, or WARmin`,
		`rep is which replacement procedure which replacement procedure is being used - either RepPlayed, RepAllAsc, RepAllDesc, RepAsc, or RepDesc`,
     `outputs a list of length 2, where the first entry says the winner(s) of the game and the second entry is how long the game took if it terminated or the cycle found if it does not terminate.`):

elif procedure=SimpPerDecks then
     print(`inputs a positive integer n and a replacement function rep`):
     print(`Finds all cycles in 2 player games with 2n unique cards where the replacement function is RepAllAsc, RepAllDesc, RepAsc, or RepDesc`):
     print(`outputs a list of all of the unique cycles found.`):

elif procedure=SimpGameLengths then
     print(`inputs a positive integer n and a replacement function rep`):
     print(`Finds the game lengths of all terminating games and the number of games that lead to cycles for all 2 player games with 2n unique cards where the replacement function is RepAllAsc, RepAllDesc, RepAsc, or RepDesc`):
     print(`Outputs a list of length 2 where the first entry is the generating function for duration of terminating games and the second entry is the number of games that lead to a cycle.`):

elif procedure=SimpGameInfo then
   print(`inputs a positive integer n and a replacement function rep`):
   print(`Finds the game lengths of all terminating games, the number of games that lead to cycles, and the all of the cycles found for all 2 player games with 2n unique cards where the replacement function is RepAllAsc, RepAllDesc, RepAsc, or RepDesc`):
   print(`Outputs a list of [Length,cycles] where Length is a list s.t. Length[1] is the generating function for duration of terminating games, Length[2] is the number of games that lead to a cycle, and cycle is a list of the unique cycles.`):

 fi:
end:



	####1. SETTING UP A GAME OF WAR###

#RanDeck(n,k,r) inputs positive integers n,k, and r. Takes a random permutation of a deck ontaining k identical copies of cards labeled {1,..,n} and simulates
#dealing out as many cards as possible such that each player has the same number of cards.
#Outputs a list of length r, where the i-th entry is a list of length floor(n*k/r) that represents the i-th players deck.

RanDeck:=proc(n,k,r) local HandSize, i, pi, Deck, decks, discard:
  Deck:=[seq(seq(i,j=1..k),i=1..n)]:
  Deck:=randperm(Deck):
  HandSize:=floor(n*k/r):
  decks:=[seq([seq(Deck[i],i=(j-1)*HandSize+1..j*HandSize)],j=1..r)]:
 end:


	###1.A. WAR PROCEDURES###
#In a war, players place $m$ cards face down and then play their next card. The following three procedures define different ways to complete a war when 
#a player does not have enough cards. 


#WARm(m,Q,r) inputs a non-negative integer m, a list Q, and a positive integer r, where Q[i] represents the i-th player's deck and |Q|=r.
#Simulates a WAR where m cards must be placed face down NO MATTER WHAT. Any player that does not have enough cards automatically loses.
#If no one can complete the WAR it outputs {0}, indicating that the game will have no winners.`):
#Otherwise, it outputs a list [tempQ,b], where tempQ is a list of each players deck AFTER the war and b is a list of the cards that were either placed face down or
#played during the WAR/.

WARm:=proc(m,Q,r) local HandSize, tempQ,i,output,a,b:
  HandSize:=[seq(nops(Q[i]), i=1..r)]:
  if max(HandSize)>m then tempQ:=Q: a:=Q:
  for i from 1 to r do 
    if HandSize[i]<m+1 then tempQ[i]:=[]:  a[i]:=[op(1..-1,Q[i]),seq(0,j=HandSize[i]..m)]: 
    else tempQ[i]:=[op(m+1..-1,Q[i])]: a[i]:=[op(1..m,Q[i])]: 
    fi: od: 
  b:=remove(t->t=0,[seq(seq(a[i][j],i=1..r),j=1..m)]):
  output:=[tempQ,b]:
  else output:={0}: fi:
  output:
 end:





WARmin:=proc(m,Q,r) local HandSize, hs, tempQ,i,output:
  HandSize:=[seq(nops(Q[i]),i=1..r)]:
  if add(HandSize[i],i=1..r)<1 then output:=0:
  else hs:=min(remove(t->t=0,HandSize)):
    if hs>m then output:=WARm(m,Q,r):
    else output:=WARm(hs-1,Q,r):
    fi:
  fi:  
  output:
 end:




WARmax:=proc(m,Q,r) local HandSize,HS, tempQ,i,output:
  HandSize:=[seq(nops(Q[i]),i=1..r)]:
  HS:=max(HandSize): 
  if HS>m then output:=WARm(m,Q,r):
   elif HS=0 then output:=0:
   else output:=WARm(HS-1,Q,r):
   fi:
  output:
end: 





	###1.B. REPLACEMENT PROCEDURES###


RepPlayed:=proc(Q,r,b,N) local c:
  c:=remove(t->t=0,b):
  [seq(Q[i],i=1..N-1),[op(Q[N]),op(c)],seq(Q[i],i=N+1..r)]:
end:

RepAllAsc:=proc(Q,r,b,N) local c:
  c:=remove(t->t=0,b):
  c:=sort(c):
  [seq(Q[i],i=1..N-1),[op(Q[N]),op(c)],seq(Q[i],i=N+1..r)]:
end:


RepAllDesc:=proc(Q,r,b,N) local c:
  c:=remove(t->t=0,b):
  c:=sort(b,`>`):
  [seq(Q[i],i=1..N-1),[op(Q[N]),op(c)],seq(Q[i],i=N+1..r)]:
end:



RepAsc:=proc(Q,r,b,N) local n,c:
  n:=nops(b)/r:
  c:=[seq(op(sort([op(i*r+1..i*r+r,b)])),i=0..n-1)]:
  c:=remove(t->t=0,c):
  [seq(Q[i],i=1..N-1),[op(Q[N]),op(c)],seq(Q[i],i=N+1..r)]:
end:


RepDesc:=proc(Q,r,b,N) local n,c:
  n:=nops(b)/r:
  c:=[seq(op(sort([op(i*r+1..i*r+r,b)],`>`)),i=0..n-1)]:
  c:=remove(t->t=0,c):
  [seq(Q[i],i=1..N-1),[op(Q[N]),op(c)],seq(Q[i],i=N+1..r)]:
end:





	####2. MULTIPLAYER SIMULATION####

OneMove:=proc(P,c,m,WAR,rep) local r, a, b, tempP, i, h, w, Q, output, W:
  r:=nops(P):
  a:=[]:
  b:=[op(c)]:
  tempP:=P:
  for i from 1 to r do
    if P[i]=[] then tempP[i]:=[0]: fi: od:
  a:=[seq(tempP[i][1],i=1..r)]:
  Q:=[seq([op(2..-1,tempP[j])],j=1..r)]:
  h:=max(op(a)):
  w:={SearchAll(h,a)}:
  if nops(w)=1 then
    b:=[op(b),op(a)]:
    output:=rep(Q,r,b,op(w)): 
  else
    W:=WAR(m,Q,r):
    if type(W,list)=true then
      Q:=W[1]:
      b:=[op(b),op(a),op(W[2])]: 
     output:=OneMove(Q,b,m, WAR,rep):
     elif type(W,integer)=true then
       output:=tie:
     else
       output:=W:
     fi: fi:
  output:
  end:



GenOneGameK:=proc(n,k,r,m, WAR,rep,K) local P, HandSize, i, Q, w:
  P:=RanDeck(n,k,r):
  HandSize:=[seq(nopsP[i],i=1..r)]:

  for i from 1 to K while nops(remove(t->t=0, HandSize))>1 and type(P,list)=true do
    Q:=P:
    HandSize:=[seq(nops(P[i]),i=1..r)]:
    P:=OneMove(P,[],m,WAR,rep):
   od:

   if i=K+1 then
    RETURN(FAIL):

   elif nops(remove(t->t=0, HandSize))=1 then
     w:=max[index](HandSize):

   elif type(P,list)=false then
     w:=P

   else
     RETURN(FAIL):
   fi:

[w,i-1]:
end:


GenManyGamesK:=proc(n,k,r,m,WAR,rep,K,w,t,M) local f,h,g,i,G,j,nw:
  for i from 1 to r do
    f[i]:=0:
    g[i]:=0:
   od:
  h:=0:
  nw:=0:

  for i from 1 to M do
    G:=GenOneGameK(n,k,r,m, WAR,rep,K):

  if G=FAIL then
    h:=h+1
   elif type(G[1],integer)=true then
     f[G[1]]:=f[G[1]]+w^G[2]
  else

   if nops(G[1])>1 then
     for j in G[1] do
        g[j]:=g[j]+t^G[2]:
       od:
    else
     nw:=nw+t^G[2]
    fi:
  fi:
od:

[[seq(f[i],i=1..r)],[seq(g[i],i=1..r)],nw,h]:
end:


GenManyGamesKEDIT:=proc(n,k,r,m,WAR,rep,K,t,M) local a, f,h,g,i,G,j,nw:
  a:=0:
  f:=0:
  g:=0:
  h:=0:
  nw:=0:

  for i from 1 to M do
    G:=GenOneGame(n,k,r,m, WAR,rep,K):

  if G=FAIL then
    h:=h+1:
   elif type(G[1],integer)=true then
     a:=a+t^G[2]:
     f:=f+t^G[2]:
  else
      a:=a+t^G[2]:
      nw:=nw+t^G[2]:
    fi:
od:

[a,f,nw,h]:
end:



#FinishGame(H,r,m, WAR,rep) inputs a list H, positive integer r where H has lenght r, non-negative integer m, a war procedure WAR, and a replacement procedure rep.
#It simulates finishing the game with r players where H[i] represents the i-th players deck and the default number of cards placed face down in a war is m,
#using WAR and rep to define how wars are played out and how cards are placed back into the winners deck.
#It terminates when the number of player's with cards left is less than or equal to 1 or a periodic orbit has been found.
#Outputs a list of length 2, where the first entry says the winner(s) of the game and the second entry is how long the game took if it terminated or the cycle found 
#if it does not terminate.


FinishGame:=proc(H,r,m, WAR,rep) local Q, P, HandSize, moves, cyc, path, N, cycle, lost, output, w:
  P:=H:
  HandSize:=[seq(nops(P[i]),i=1..r)]: 
  moves:=0:
  cyc:=0:
  path:=[P]: 

  while nops(remove(t->t=0, HandSize))>1 and type(P,list)=true do
    P:=OneMove(P,[],m,WAR,rep):
    HandSize:=[seq(nops(P[i]),i=1..r)]:
    moves:=moves+1: 

    if moves mod 10=0 then
      N:=Search(P,path): 

      if N>0 then
        cyc:=nops({op(N..moves,path)}):
        cycle:=[op(N..N+cyc-1,path)]:
        
        HandSize:=[seq(nops(P[i]),i=1..r)]:
        lost:={SearchAll(0,HandSize)}:
        w:={seq(i,i=1..r)} minus lost:
       break:
     fi: fi: 
    path:=[op(path),P]:
    od:

    if cyc>0 then
      output:=[w,cycle]:

     elif nops(remove(t->t=0, HandSize))=1 then
       output:=[max[index](HandSize),moves]:

     elif type(P,list)=false then
       output:=[P,moves]:

     else
       RETURN(FAIL):
     fi: 

  output:
end:



	####3. FINDING CYCLES IN SIMPLE WAR GAME####

#SimpPerDecks(n,rep) inputs a positive integer n and a replacement procedure rep. 
#Outputs a list of unique cycles found in two player games with 2n unique cards and replacement procedure rep.

SimpPerDecks:=proc(n,rep) local Decks, i, counter, LedToCyc, cycles,o,i1,j1,j2, C,N, H,G,a,b,k,M,N1, DrZ:
 if rep<>RepAsc and rep<>RepDesc and rep<>RepAllAsc and rep<>RepAllDesc then
    print(`This is not a simple war game. Can't use this procedure.`):
  else
    #Here, we do not need to distinguish between Player 1 and Player 2 to figure out how all possible games will play out.
    #Because the order that cards are replaced is the same for Player 1 and Player 2, running the simulation where player 1 starts of with the deck H[1] and 
    #Player 2 starts off with the deck H[2] will tell us everything we need to know about the game where player 1 starts with H[2] and player 2 starts with H[1].
    #Because the player with the (unique) highest card will never lose that card, there will never be a point in the game where player 1 and player 2 swap decks.
    #So, running simulations where we look at the partitions of 2*n into two parts of size n where player 1 always starts with the highest card, we get also get half 
    #of the unique cycles (to obtain the other half, just switch player 1 and player 2's decks).
 
    Decks:=setpartition([seq(i,i=1..2*n)],n):  
    #It is obvious that if one player's highest card is lower than the other player's lowest card, then they lose the game in n rounds.
    #So, we do not need to run any simulations for the first partition in Decks.


    #Right now, Player 1 always starts with 1. For later in the code, it is easier if he starts with 2*n since he can never lose this card.
    #We need to turn all of the current 2*n's into 1's and all of the current 1's into 2*n's. Note, that right now the current 1's are always 
    #the first element of P1's deck
    Decks:=subs(2*n=DrZ,Decks):  
    Decks:=subs(1=2*n,Decks):
    Decks:=subs(DrZ=1,Decks): 
    
    #It is obvious that if one player's highest card is lower than the other player's lowest card, then they lose the game in n rounds.
    Decks:=remove(t->t=[[2*n, seq(i,i=2*n-n+1..2*n-1)],[seq(i,i=2..2*n-n),1]],Decks):

    N1:=n!: 
    N:=nops(Decks): 
    counter:=0:
    LedToCyc:=0:
    cycles:=[]:
    C:=[seq([],i=1..2*n)]: 
    o:=permute(n):
   
    #Right now each players deck is in a specific order, so we need to look at the various ways to shuffle them.
   for i1 from 1 to N do
    for j1 from 1 to N1 do
      for j2 from 1 to N1 do
        H:=[Decks[i1][1][o[j1]],Decks[i1][2][o[j2]]]: 
        
      if counter mod 50000 = 0 and counter>0 then     	   
	print(`Played`, counter, `games`);
	print(nops(cycles), `distinct cycles`);
	print(LedToCyc, `got a cycle`);
       fi:

      counter:=counter+1:

      G:=FinishGame(H,2,1, WARm,rep):
      #Note: m and WAR don't matter since every card is unique.

      if type(G[2],list)=true then 
         #We have a cycle!
         LedToCyc:=LedToCyc+1:
         a:=G[2]:  
         b:=[op({op(a)})][1]: 
         k:=b[1][2]:  
         M:=Search(b,C[k]):

        if M=0 then 
          C:=[seq(C[i],i=1..k-1), [op(C[k]),b], seq(C[i],i=k+1..2*n)]: 
          cycles:=[op(cycles),a]:
         fi: fi:
      od: od: od:
    print(2*nops(cycles), `distinct cycles`):
    print(2*LedToCyc, `got a cycle`):
    print((2*n)!, `total games`):
  fi:
	#We need to consider when Player 1 and Player 2 switch decks.
       [op(cycles), seq(cycles[i][2],cycles[i][1],i=1..-1)]:
  end:




#SimpGameLengths(n,rep) inputs a positive integer n and a replacement procedure rep. 
#Outputs a list of length 2 where 
#the first entry is the generating function for duration of terminating 2 player games with 2n unique cards and replacement procedure rep
#and the second entry is the number of games that lead to a cycle.

SimpGameLengths:=proc(n,rep) local Decks, i,o,i1,j1,j2, counter, LedToCyc,N, f, H,G,N1,DrZ:
 if rep<>RepAsc and rep<>RepDesc and rep<>RepAllAsc and rep<>RepAllDesc then
    print(`This is not a simple war game. Can't use this procedure.`):
  else
    #Here, we do not need to distinguish between Player 1 and Player 2 to figure out how all possible games will play out.
    #Because the order that cards are replaced is the same for Player 1 and Player 2, running the simulation where player 1 starts of with the deck H[1] and 
    #Player 2 starts off with the deck H[2] will tell us everything we need to know about the game where player 1 starts with H[2] and player 2 starts with H[1].
    #Because the player with the (unique) highest card will never lose that card, there will never be a point in the game where player 1 and player 2 swap decks.
    #So, running simulations where we look at the partitions of 2*n into two parts of size n where player 1 always starts with the highest card, we get also get half 
    #of the unique cycles (to obtain the other half, just switch player 1 and player 2's decks).
 
    Decks:=setpartition([seq(i,i=1..2*n)],n):  


    #Right now, Player 1 always starts with 1. It is easier if he starts with 2*n since he can never lose this card.
    #We need to turn all of the current 2*n's into 1's and all of the current 1's into 2*n's. Note, that right now the current 1's are always 
    #the first element of P1's deck
    Decks:=subs(2*n=DrZ,Decks):  
    Decks:=subs(1=2*n,Decks):
    Decks:=subs(DrZ=1,Decks):
    
    #It is obvious that if one player's highest card is lower than the other player's lowest card, then they lose the game in n rounds.
    Decks:=remove(t->t=[[2*n, seq(i,i=2*n-n+1..2*n-1)],[seq(i,i=2..2*n-n),1]],Decks):

    N1:=n!: 
    N:=nops(Decks): 
    counter:=0:
    LedToCyc:=0:
    f:=0:
    o:=permute(n):
   
    #Right now each players deck is in a specific order, so we need to look at the various ways to shuffle them.
   for i1 from 1 to N do
    for j1 from 1 to N1 do
      for j2 from 1 to N1 do
        H:=[Decks[i1][1][o[j1]],Decks[i1][2][o[j2]]]: 
      if counter mod 50000 = 0 and counter>0 then     	   
	print(`Played`, counter, `games`);
	print(nops(cycles), `distinct cycles`);
	print(LedToCyc, `got a cycle`);
  	print(f):
       fi:

      counter:=counter+1:

      G:=FinishGame(H,2,1, WARm,rep):
      #Note: m and WAR don't matter since every card is unique.

      if type(G[2],list)=true then 
         #We have a cycle!
         LedToCyc:=LedToCyc+1:

      elif type(G[2],integer)=true then
         f:=f+t^G[2]:
        fi:
      od: od: od:
    print(2*LedToCyc, `got a cycle`):
    print((2*n)!, `total games`):
  fi:
  #We need to consider when Player 1 and Player 2 switch decks, as well as the 2(n!) games where one player's highest card is lower than the other player's lowest card.
      [2*f+2*n!^2*t^n,2*LedToCyc]:
  end:

        
#SimpGameInfo(n,rep) inputs a positive integer n and a replacement procedure rep. 
#Outputs a list of length 2 where 
#the first entry is a list of length two whose 
	#first entry is generating function for duration of terminating 2 player games with 2n unique cards and replacement procedure rep
	#and the second entry is the number of games that lead to a cycle.
#and the second entry is the list of unique cycles.

SimpGameInfo:=proc(n,rep) local Decks, i,o,i1,j1,j2, counter, LedToCyc, cycles, C,N, f, H,G,a,b,k,M,N1,DrZ:
 if rep<>RepAsc and rep<>RepDesc and rep<>RepAllAsc and rep<>RepAllDesc then
    print(`This is not a simple war game. Can't use this procedure.`):
  else
    #Here, we do not need to distinguish between Player 1 and Player 2 to figure out how all possible games will play out.
    #Because the order that cards are replaced is the same for Player 1 and Player 2, running the simulation where player 1 starts of with the deck H[1] and 
    #Player 2 starts off with the deck H[2] will tell us everything we need to know about the game where player 1 starts with H[2] and player 2 starts with H[1].
    #Because the player with the (unique) highest card will never lose that card, there will never be a point in the game where player 1 and player 2 swap decks.
    #So, running simulations where we look at the partitions of 2*n into two parts of size n where player 1 always starts with the highest card, we get also get half 
    #of the unique cycles (to obtain the other half, just switch player 1 and player 2's decks).
 
    Decks:=setpartition([seq(i,i=1..2*n)],n):  
    #It is obvious that if one player's highest card is lower than the other player's lowest card, then they lose the game in n rounds.
    #So, we do not need to run any simulations for the first partition in Decks.


    #Right now, Player 1 always starts with 1. For later in the code, it is easier if he starts with 2*n since he can never lose this card.
    #We need to turn all of the current 2*n's into 1's and all of the current 1's into 2*n's. Note, that right now the current 1's are always 
    #the first element of P1's deck
    Decks:=subs(2*n=DrZ,Decks):  
    Decks:=subs(1=2*n,Decks):
    Decks:=subs(DrZ=1,Decks): 
    
    #It is obvious that if one player's highest card is lower than the other player's lowest card, then they lose the game in n rounds.
    Decks:=remove(t->t=[[2*n, seq(i,i=2*n-n+1..2*n-1)],[seq(i,i=2..2*n-n),1]],Decks):

    N1:=n!: 
    N:=nops(Decks): 
    counter:=0:
    LedToCyc:=0:
    cycles:=[]:
    C:=[seq([],i=1..2*n)]: 
    f:=0:
    o:=permute(n):
   
    #Right now each players deck is in a specific order, so we need to look at the various ways to shuffle them.
   for i1 from 1 to N do
    for j1 from 1 to N1 do
      for j2 from 1 to N1 do
        H:=[Decks[i1][1][o[j1]],Decks[i1][2][o[j2]]]: 
        
      if counter mod 50000 = 0 and counter>0 then     	   
	print(`Played`, counter, `games`);
	print(nops(cycles), `distinct cycles`);
	print(LedToCyc, `got a cycle`);
  	print(f):
       fi:

      counter:=counter+1:

      G:=FinishGame(H,2,1, WARm,rep):
      #Note: m and WAR don't matter since every card is unique.

      if type(G[2],list)=true then 
         #We have a cycle!
         LedToCyc:=LedToCyc+1:
         a:=G[2]:  
         b:=[op({op(a)})][1]: 
         k:=b[1][2]:  
         M:=Search(b,C[k]):

        if M=0 then 
          C:=[seq(C[i],i=1..k-1), [op(C[k]),b], seq(C[i],i=k+1..2*n)]: 
          cycles:=[op(cycles),a]:
         fi: 

      elif type(G[2],integer)=true then
         f:=f+t^G[2]:
        fi:
      od: od: od:
    print(2*nops(cycles), `distinct cycles`):
    print(2*LedToCyc, `got a cycle`):
    print((2*n)!, `total games`):
  fi:

  #We need to consider when Player 1 and Player 2 switch decks, as well as the 2(n!) games where one player's highest card is lower than the other player's lowest card.
      [[2*f+2*n!^2*t^n,2*LedToCyc],[op(cycles),seq(cycles[i][2],cycles[i][1],i=1..-1)]]:
  end: